Practical wavelet design on the sphere 
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Abstract 

We address the question of designing isotropic analysis functions on the sphere which 
are perfectly limited in the spectral domain and optimally localized in the spatial 
domain. This work is motivated by the need of localized analysis tools in domains 
where the data is lying on the sphere, e.g. the science of the Cosmic Microwave 
Background. Our construction is derived from the localized frames introduced by 
Narcowich et al. (2006). The analysis frames are optimized for given applications 
and compared numerically using various criteria. 



Introduction 



Localized analysis for spherical data has motivated many researches during 
the past deca de. Data defined on the sphere are studied in domains as variou s 
as cosmolo gy (Hinshaw et al. . 20061: [liivon et aL . 20021: McEwen et al. . 2007 ). 
geophysics ( Holschneider et al. . 20031 : Wieczorek and Simon I l200nh . medicine, 
computer vision. When dealing with data on the whole sphere, spectral analy- 
sis can be achieved by Spherical Harmonics Transform (SHT) - the equivalent 
of the Fourier Series on the circle. But in many practical situations, data are 
defined or available on a subset of the sphere only. For example, cosmologists 
try to give sharp estimates of the cosmic microwave background (CMB) or 
its power spectrum but strong foreground emissions superimpose to the CMB 
making the observations unreliable for CMB studies. Moreover, fully observed 
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clean non stationary fields or stationary fields with additive non-stationary 
noise still require spatially localized tools. In such situations, the SHT is not 
adequate, because of the poor spatial localization of the basis functions. In the 
case of Euclidean spaces, in which the Fourier Transform suffer from the same 
lack of localization, multiscale and wavelets theory provide a mathematically 
elegant solution of proven practical efficiency. 



Adaptation to the sphere of the "wavelet" transform (in the broad sense of fil- 
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). Since then, Antoine & Vandergheynst (1999) showed that any Continu- 



ous Wavelet Transform (CWT) on the sphere can be viewed locally as a regular 
CWT on the Euclidean tangent planes, th anks to the stereographic corr e spon- 



dence between the sphere and the plane (lAntoine and VandergheynstI , Il999l : 



Wiaux et al.l . 120051 ). One can then "forget" the sphere by projecting it on tan- 
gent planes, realizing the analysis in these planes, and then apply the inverse 
projection to get back eventually to the sphere. A discreti z ed ve rsion of this 



approach of CWT has been presented by Bogdanova et al. (l2005l ). leading to 



wavelet frames. This approach has already been followed in astrophys i cs for 



the analysis of t he Co smic Microwave Background (CMB) (IVielva et al.l . 12004 : 



McEwen et al.l . 120071 ) . However these wavelets are usually defined in the spa- 
tial domain and have infinite support in the frequency domain (which must 
be truncated in practice). 



In the present work, we follow and extend the approach of iNarcowich et al. 
(I2OO6I ) and their construction of "needlets". A similar construction can be 



found in IStarck et al.l (|2006l ). The needlet transform has important charac- 
teristics. Firstly it is intrinsically spherical. No intermediate tangent plane is 
needed to define it. Secondly, it does not depend on the particular spherical 
pixelization chosen to describe the data. Thirdly, although the needlets still 
have an excellent spatial localization, they have a finite spectral support ad- 
justable at will . They are axisymmetric (which is convenient when dealing 
with statistically isotropic random fields) and thus the needlet coefficients are 
easily computed in the Spherical Harmonics (Fourier) domain. Data filtering is 
defined by multiplication of the Spherical Harmonics coefficients by well cho- 
sen window functions (which is equivalent to convolution in spatial domain). 
Needlets are well defined in theory and the statistical properties of their coeffi- 
cients have already been established for isotropic Gaussian fields f Baldi et al.l 
(|2005). However, the performance of a needlet-based analysis depends on the 
particular shape of the needlet. 



This paper considers spherical filters which are generalizations of needlets in 
the spirit of dual (non-tight) analysis and reconstruction frames. We focus on 
the design issue, namely the optimization of the window functions (that define 
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the isotropic filtering operations) for some given tasks. We consider only band- 
limited needlets. This choice is motivated by applications in high-precision 
cosmology. Indeed, the CMB power spectrum is highly dynamic (few peaks and 
power-law decay) and good subsequent cosmological parameters estimation 
requires high accuracy in some critically delimited spectral ranges. Once the 
range is fixed, we optimize the shape of window functions in two directions: 
1) By requesting the best spatial localization of associated needlets, in an 
en ergv-sense (L^) whic h is easily solved. This is an application of the work 
of Simons et al. (I2OO6I ) which adapted to the sphere the problem solved by 



Slepian " (jl978l ) on the real line, giving rise to the well known prolate spheroidal 
wave functions (PSWF). 2) By following statistical considerations: given some 
region ("mask") in which the data is missing or thrown away and assuming that 
the full data is the realization of some Gaussian isotropic random field (this is 
the usual assumption made on the CMB), we minimize the mean integrated 
square error due to the mask in the needlet analysis outside the mask. More 
criteria and applications to cosmological science will be given in a future work. 

The paper is organised as follows. In Section [H we expose the general con- 
struction of needlets. In Section [21 we define and optimize the two criteria 
(geometrical and statistical) which provide localized analysis filters. Their ef- 
ficiency is illustrated in Section [3] with numerical simulations following the 
model of a masked observation of the CMB. The proofs are postponed to 
Appendix [Al 



1 Needlets frames 



1.1 Background and notations 



Denote § the unit sphere in with generic element ^ = {6, ip) in spherical 
polar coordinates: 9 G [0, tt] is the colatitude and ip G [0,27r[ the longitude. 
Let H = L^(S) be the space of complex- valued square integrable functions 
on S under the Lebesgue measure d^ = sin6'd6'dv9. Endowed with the inner 
product {f,g) := Is f{09*i0^^i H is a Hilbert space. Let || ■ || denote the 
associated norm on M. The usual complex spherical harmonics on S (which 
definition is recalled in Appendix [B1) are denoted (y£m)^>o,-^<m<£- They form 
an orthonormal basis of H. 

In the following, we consider a field X G H. Its random spherical harmon- 
ics coefficients or multipole moments are denoted agm = {X, Ygm) ■ H can be 

± 

decomposed in harmonic subspaces: H = H^, where is the linear span 
of Yim,m = —£, ■ ■ ■ ,i. The number i is referred to as the multipole number 
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or frequency (understood as a spatial frequency). Let Hi be the orthogonal 
projection on Mi. It has an expression involving Legendre polynomials Li (see 
Appendix [B| 

UiXio = E {X, Yi^)Yim.{0 = J xiOHi ■ e')de'- (1) 

where = cos 9 cos 9' + sin 9 sin 9' cos{(p — ip') is the usual dot product on 



A mapping on S which depends on the colatitude 9 only is said to be axisym- 
metric. The convolution of a bounded axisymmetric function H{^) = h{cos9) 
with an arbitrary spherical function X is well defined through 

H*x{o= f (2) 

The convolution theorem holds: 

H *X = Y^ hiaimXim- (3) 

Im 

where aim = Ytm) are the multipole moments of X and hi are the Legendre 
series coefficients of i.e.h = J^em^iLi. Then, an isotropic wavelet analysis 
can be implemented either in the spatial (i.e. direct) domain using ([2]) or in the 
harmonic domain using (JSj). We choose the latter, which accounts to multiply 
the harmonic coefficients of the field of interest X by a spectral window (hi). 
For a countable index set JT, let {h^^^)j^j be a family of window functions 
in harmonic domain : h^^^ E £°°(N). The corresponding harmonic smoothing 
operators on HI are 

vl>(^) = Y: h?Ui. (4) 

We call exact reconstruction condition the one ensuring that E = Id. It 
also writes 

E h^'^ = 1 (5) 

In the following, j is referred to as the scale, in analogy with the multiresolu- 
tion analysis terminology. Important examples of windows families having the 
property (l5|) may be obtained thanks to the i?-adic mechanism: let S > 1, 
i7 = { — 1}UN, h^"^^ = 5Q{i) and the spectral windows be all related to a 
continuous function h by 

Vj G N, = h ( Jy) • (6) 

If h is compactly supported on [-^, B], then each window h^^^ may overlap with 
adjacent windows h^^~^^ and h'3+^'> only. The exact reconstruction condition 
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in this case is satisfied as soon as 

Vx e [1,5], h(x) + h{B-^x) = 1 
This example is illustrated in Figures [T] and [2l 



(7) 





(a) Original map 




(b) Smoothed maps, scales j = 2, 5 




(c) Smoothed maps, scales j = 6, ...,9 

Fig. 2. Input map of a CMB sky (from WMAP), and corresponding smoothed maps 
(with the spline filters of Figured]). 
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1.2 Needlet tight frames 



Recall that a countable family of functions {/„} in a Hilbert space 7i is a 
frame with frame bounds Ci, C2 if 



\/9 en , C,\\g\\H<Y.\i9Jn)n\' <C2\\g\ 



It is a tight frame if we can choose Ci = C2. Frames can be thought of redun- 
dant "bases", and this redundancy can be exploited for robustnes s issues. The 
tight ness property is valuable in terms of numerical stability (see IPaubechied . 
19921 . Chap. 3 and the references therein). 



The constructi on that follows is fr om lNarcowich et al.l ( 120061 ). The term needlet 
was coined by iBaldi et al.l (120061 ). Let K he a finite index set and {^k}keK G 
§1^1 a set of quadrature points on the sphere, associated with a set {\k}k£K £ 
M'^I of quadrature weights. 

Definition 1.1 (Quadrature) {{C,k, Xk)}keK is said to provide an exact Gauss 
quadrature formula at degree £max if 



e=o keK 
This quadrature formula is said positive-weight if Xk > 0,k G K. 



Remark 1.2 We refer to Woroshkevich et al\ l200a ) for an example of a 
proper choice of quadrature points and weights that fulfils this property ( called 
GLESP). Other pixelization schemes such as HEALPix jiGdrski et all . \200A) 
fulfil approximately this property with a number of points of order Cll^^^ and 
quadrature weights of order jrk — /or some positive constant C. 



Suppose that the window functions h^^^ are non-negative and with finite spec- 
tral support. Define 

e N, b^p := (8) 

and d^^^ := max{£ : h^/^ ^ 0} (in the 5-adic case, d^^^ = B^^^). For each scale 
j, we have a pixellization {ik\x\^^}f^^j^(j). 

Definition 1.3 (Needlets and Needlet coefficients) For every j G J and 
every index k G K'^^'^ the function 



(9) 
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is called a needlet. For X G M, the inner products {X^ip^^^) are called needlet 
coefficients and are denoted I3^\ 

Up to a rotation of the sphere putting ^^■'^ on the North pole and to the 
multiplicative term \l A^f •* , all the needlets of a given scale j have exactly the 
same shape. In particular, they are axisymmetric. When (. ^ is sufficiently 
smooth, one gets the intuition from ^ that the needlet -^ll^^ is localized around 

if- 

The following Proposition state that the harmonic smoothing operation de- 
fined by (111) can be seen as the decomposition of H on the needlets family 
built with (IHl) , and that this fa mily is a tight frame. It is a straightforward 



adaptation of iBaldi et al.l (120061 . Proposition 2.3) 



Proposition 1.4 Let j e J . Assume that {{^k \ ^k^)}keK(i) provides an ex- 
act and positive-weight quadrature formula at degree 2d^^K Then 

ksKU) 

Assume that for any j G J^, {{^k\ ^k^)}keKO) provides an exact and positive- 
weight quadrature formula at degree 2d^^\ Under the exact reconstruction con- 
dition 

VXGH, X=E E and ||Xf = ^ ^ l^^'^^ • 



Remark on Terminology The analysis of an input field X in the way 
described above is called filtering. This filtering has two equivalent expressions, 
in the spatial and in the spectral domains; see the convolution formula ([3]). 
These expressions involves two "dual" mathematical objects : the functions 
h''^^ and of the frequency called window functions or spectral windows, 
and the spherical functions if)^^^ called needlets, which are nothing else but 
the rotated axisymmetric functions built from the Legendre transform of h^^^ 
(see Definition II. Sp . We call filter either of the two above objects, when the 
domain (spatial or spectral) is not specified. 



1.3 Generalized needlet frames 



We are concerned with the development of a flexible spectral analysis on the 
sphere which remains practical at high resolution. The forecoming CMB ex- 
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periment PlancJ^ will provide 50 mega-pixel maps with accuracy such that 
multipole moments will be reliable up to £ ^ 4000. 

For maximum flexibility, we shall consider constructions which are not neces- 
sarily dyadic nor i?-adic. This is motivated by applications, as described in the 
Introduction. Moreover, we will design analysis frames which will not be neces- 
sarily tight. Their dual frames will be the corresponding reconstruction frames. 
This allows fine tuning of the localization properties of the decomposition func- 
tions but it is also well known that it does not ensure similar properties for 
the reconstruction functions. Nevertheless, for the application goals discussed 
in the introduction, we will design strictly band-limited needlets with support 
:= [^HL.^^^L]. ^min > if j > 0. Then the subsequent "wavelet design" 
operations will be performed in the harmonic domain. 

Since the needlet coefficients jS^^^ and P^V of a Gaussian stationary {i.e. 
isotropic) field are independent if L^^^ fl L^^ ^ = 0, the bands L^^^ are cho- 
sen to overlap as little as possible. Other choices are possible; for instance 



Starck et al.l (120061 ) take overlapping spectral windows supported on [0, 2-']. 



The three ingredients for our spherical "multi-resolution" approach are harmonic- 
space implementation, dual wavelet frames and spectral window design. In this 
subsection, we briefiy describe the first two elements. In Section [2, we go into 
the theory and practice of window design. 



1.3.1 Dual frames 

Proposition II .41 shows that the needlets of Definition 11.31 with ([8]) can be used 
in both analysis (or decomposition) and synthesis (or reconstruction). This 
accounts to say that the needlet frame is its own dual frame. We choose to 
keep the Definition 11.31 of the needlets and associated coefficients but to relax 
condition ([8]). By sacrificing the tightness of the frame, we gain much freedom 
in the design of the spectral windows. Also, the precise space-frequency picture 
provided by the needlet construction is preserved. 

From any windows family {W^)j(zj such that V£ G N, Z^jej i^^P^ > Oi define 
the synthesis windows by 

Vj G J, V£ G N, hP = (10) 

and put h^^^ := W^h^^'^ so that ([5]) easily follows. We retain Definition 11.31 
for the decomposition needlets and needlets coefficients and further define the 

^ see www.rssd . esa. int/Planck/. 
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reconstruction needlets as 

^F(0 = \/AFE&?Me-e?). (11) 

1=0 

Proposition 1.5 Assume that there exists positive constants Ci, C2 such that 

V£ G N,Ci < E < C'2 . (12) 

Assume that for any j G J, the set {{ik\'^^k'')}k&Kij) provides an exact and 
positive-weight quadrature formula at degree 2d^^\ Then, under the exact re- 
construction condition the family {ip^^} is a frame with frame bounds 
constant Ci and C2. Its dual frame is the family {ipk^}- In particular 

VX G H, X i 5: (3m^ and ||Xf = J] ^ ^l^'^(3l^'> , (13) 

withp^i^ ■.= {x,^^i^). 

Define the analysis, synthesis and smootliing operators at scale j G ^7 by 

= Y.ihflli, I'(J) =Y,hflli and ^(^^ = <|(j)<l>(j), respectively. Then, the 

i 

exact reconstruction formula ^''"'^ = M holds true. 

An example of an analysis/synthesis windows family following this scheme is 
displayed in Figure [3l in which we took optimally concentrated PSWF (see 
Section [21) functions for analysis. It illustrates the fact that this choice does 
not lead to well localized synthesis needlets (as their spectral shapes are non 
smooth). However, this may not be a shortcoming if one is interested in the 
needlet coefficients (3^"^ = (X, tp^^^) per se, which reflect the local properties of 
the fleld X. 



1.3.2 Practical computation of needlet coefficients 

Evaluation of inner products (X, ipl^^) in the direct space is practically unfeasi- 
ble from a pixelized sphere at high resolutions. The needlet coefficients (3^^^ are 
thus computed via direct and inverse harmonic transforms as a consequence 
of the following Proposition. 

Proposition 1.6 The needlet coefficients verify 4^'^ = \/AF'^^^^^(d^^)- 

The computation of the smoothed field ^'^•'^X is performed in the harmonic 
domain by multiplying the multipole coefficients a^m of X by the factors b^^^ . 
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Multipole 



Fig. 3. i?-adic analysis (top) and corresponding synthesis (bottom) window functions 
(j = 6,...,ll ; B = l.7). 

Finally, the needlet coefficients are retrieved as the values of ^^^^X at the 
points ^[^^ up to a multiplicative term. Starting from the field X sampled at 
some quadrature points, this operation is summed up by the diagram 

whereas the synthesis operation is summed up by 



Standard pixelization packages, such as HEALPix, GLESP or shtoolS 
come with optimized implementations of the direct and inverse Spherical Har- 
monic Transforms. For example, in the HEALPix scheme, pixels are located on 
rings of constant latitude, allowing for fast SHT. This makes the computation 
easy and tractable even at high resolution. The needlet coefficients at a given 
scale j can be visualized as a pixelized map. If the quadrature weights {A^"'^} 
are equal, the smoothed maps of Fig. [21 which are the outputs of the processing 
( |T4l) , provide a precise and easily interpretable picture of the space- frequency 
analysis. 



available at http:/ /www. ipgp.jussieu.fr/~wieczor/SHTOOLS/SHTOOLS. html 
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Remark 1.7 The quadrature points and weights {{Q , )}k£KU'i ^"^'^ de- 
fine the needlet coefficients and to sample the smoothed field "^^^^X may be 
chosen identical to {{^k, ^k)}keK used to sample the input field X. However, 
for data compression and computational efficiency, one can consider alterna- 
tively to take the minimal K^^^ providing an exact positive-weight quadrature 
formula at a proper degree. 



2 Design of optimally localized wavelets 



In this section, we define some criteria to compare the window profiles. Some 
of them are easily optimized, others are only investigated numerically. We 
first give some examples of generic needlet profiles we can think of (Sec- 
tion EH]). Then, we restrict ourselves to a single scale j and an associated 
band L := [£min,^max]- The superscript (j) will be omitted in the notations 
when no confusion is possible. We present the (Section [2. 2p and statistical 
(Section [2.3p criteria, with practical implementation details on their optimiza- 
tions. 



2. 1 Examples. 



Narcowich et al.l (120061 ) have derived the following theoretical bound that con- 
trols the decay of the needlets. In the S-adic case, if the function b := 
defining the analysis spectral window is M-times continuously differentiable, 



1 + LBJ-i arccos(^ ■ 



for some constant C = C(b). This condition still allows a wide range of possi- 
bilities for designing the function b. Without restricting ourselves to the -B-adic 
case, we implemented solutions to optimize in practice, non asymptotically, 
the shape of windows b^^^ regarding some applications. 

To illustrate the kind of aspects we are concerned with, we compare in Figured] 
the azimuthal profiles (in the spatial domain) of various axisymmetric needlets. 
The needlets are built from window functions b^^^ via relation (J9]) and = 
(0,0), i.e. they are centered on the North pole, and then are considered as 
functions of 9 only. This illustration is restricted to the 9**^ dyadic scale, i.e. 
frequencies in the band L := [256,1024]. We shall compare heuristically five 
families of window functions. Note that the last two are not limited to band 
L. 
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(1) Square roots of splines of various orders. For any odd integer M, there 
exists a spline function h of order M, non-negative, compactly supported 
on [|,2] and such that the h\^^'s defined by (I6|) verify Q. It remains to 



define b'f^ = ^Jh'f\ 

(2) Best concentrated Slepian functions in caps of various radii (cf Sec- 
tion [22D- The window function fc^"*^ is the minimizer of the criterion (i20l) . 
It is band-limited on L and optimally concentrated in a polar cap = 
{C, : 9 < Oq}), Oq being a free parameter. 

1 

(3) Denote G a primitive of the C°° function g : x t-^ e T^l(_i^i)(a;) and 
put 

b(x) = G(-8x + 3) -G(-4x + 3) (15) 
and b\^^ = b (^) . This window function is used in Pietrobon et al. ( 2006h . 



(4) From the i?-spline function of order 3 

Bs{x) = - 2|^ - 4|x - 1|^ + 6|a;|^ - 4\x + 1|^ + \x + 2|^), (16) 

form b(x) = |(-B3(2x) — Bs{x)) and define b^p = b(^). This window 



function is used by IStarck et al.l (120061 ) 



The Mexican hat wavelet on the sphere is the function the stereographic 
projection of which on the Euclidean plane is the usual Mexican hat 
wavelet. It has the following close expression depending on some positive 
scale parameter R 

iJaiO) oc (1 - 2i?2 tan2(^/2)) exp{-2i?2 tan2(^/2)}. (17) 



This wavelet is popular in the astrophysics community (see e.q. lGonzalez-Nuevo et al. 



20061 ). We have chosen i? = 6.10 ^ such that the spectral window is al- 



most zero for £ > 1024. 



2.2 IL -concentration and variations 



Our first attempt to achieve a good spatial localization of a needlet is to 
optimize a L^-norm based criterion, adapting to the sphere a problem that 
is well-known on the real line. In their seminal work in the 1960s and 1970s, 
Slepian and his collaborators studied the properties of prol ate sphe r oidal wave 
functions (PSWFs) in the ID case of real functions (see Slepian . 1983i . and 



the references therein). PSWFs may be defined as functions with optimal 
energy concentration in the time domain, under some band-limitation con- 
straint. Equivalently, they are the eigenfunctions of a time-frequency concen- 
tration kernel or the solutions of a Sturm-Liouville differential equation. The 
time-frequency concentration of PSWFs is understood in terms of continuous 
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(b) PSWFs localized in polar caps of 0.5, 1, 1.5 and 5 degree opening 
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(c) (red curve) Exponential function described in Eq. (fT5]l . (green curve) S-spline 
function of Eq. (fT6l) and (blue curve) Mexican hat described in Eq. (fTTl) . 



Fig. 4. In left column, the shape of the spectral windows as a function of i. In 
middle and right columns, the profile of the filters is plotted in the spatial domain 
as a function of 6 {6 in degrees) with logarithmic and linear scales respectively, to 
illustrate both the decrease of the tail of the needlets far from the North pole and 
the shape of their first bounces. 

Fourier transform on M. A dis crete ve r sion o f this theory, based on Fourier 
series coefficients, is derived in ISlepianI (Il978l ). 



In the last few years, W alter and coauthors exp loited th ese ID PSWFs to 



deriv e Slepian series (in IWalter and Shenl . l2003l : see also 
2004 ) , and wavelets based on th e best concentrated PSWF (jWalter and Shenl . 
20041 : IWalter and Soleski l2005l ). 



Moore and Cada 



On the sphere, we shall on ly consider the equivalent of Discrete PSWFs, fol- 
lowing lSimons_eLalJ (120061). From a window function {hi} with support L, 
define the axisymmetric function i}) by 



(18) 



The set of functions ip of the form (|T8l) is denoted Bl C H. Given a spherical 
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domain consider the minimization, among non-zero functions in Bl, of the 
criterion 

This exte nsion to the sphere o f Slepian's concentration problem is studied in 
details by Simons et al. (I2OO6I ) in the case £min = 0. We call PSWF (by abuse 



of language) and denote i/Jq a normalized minimizer for Cq{iP). 



The criterion (fTOl) has a simplified expression when ^2 is axisymmetric. Con- 
sider the polar cap Qe^^ = '■ < 9o} and define the coupling matrix 

D = {Di^e>)i^i>(zL by 



= I Le{z)Li>{z)dz 

'4-1) "'cos6»o 



and 



Then 



Hi;) 



CnW = 1 



2^max ~l~ 1 , 



b*Db 



(20) 



and the minimization of ( fTOl) becomes an eigenvalue problem. The solution of 
this minimization depends on the opening Oq. In Figure [5] we plot the value 
of Cug^ against Oq for ip^^^ , ip^^^ , ip^^^^ . The lowest curve is the minimum of 
the criterion for all openings 6*0. It is clear that there is no optimal function 



uniformly in 9q: the concentration criterion Cqq 
the best possible value for ^0 = ^1 only. 



of each PSWF ip^ reaches 



o 
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o 
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Fig. 5. Localization for L^-energy criterion of PSWFs, band-limited into L = [33, 64]. 
The dashed line is the minimum of the criterion Cfig^ as a function of and the 
other ones are the values of Cn^^ (ip) evaluated at ip = ipQ^ , tpQ^ and ipti-^g ■ 
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As in the 1-dimensional case, the spectrum of D exhibits a "step function" 
behaviour: denoting = tr D (the "Shannon number"), the matrix D has 
about e igenvalues very clos e to 1, and most of the others close to zero (see 
FigEl and lSimons et al.ll2006l for details). 




Fig. 6. Eigenvalues of matrix D with 
number N = 13.3. 



50° and L = [17, 64]. In this case, Shannon 



When several eigenvalues of D are extremely close to 1, it is computation- 
ally difficult to find the largest one and the associated eigenvector. In the 
case of Q a polar cap and i^ain = 0, one can advantageously solve the less 
degenerat ed eigenvalue problem associated with the Griinbaum differential 



equation ( IGriinbaum et al.l . 119821 ) which has the same solutions as (ITQll . We 
are not aware of an equivalent theory in the case ^min > 0. 

With e being of the order the machine precision, all vectors in K = Ker(D— 

A>l-e 

Aid) have well spatially localized counterparts, but they are not necessarily 
positive (in harmonic domain). This is not acceptable for instance if we were 
to use them as windows associated to smoothing operator (denoted h in the 
first Section), and implement this operator using a needlet analysis-synthesis 
scheme, the window of which has to be defined as the square-root of the 
PSWF's window. To circumvent this, we therefore introduce a modified cou- 
pling matrix D = D + aH*H where a > is a tuning parameter and H is the 
tridiagonal second-order finite difference matrix. Window functions are now 
obtained as minimizers of Cq(^) = 1 — 'jj^^ instead of Cf^. The additional term 
favors non-oscillating functions among the vectors of which are undistin- 
guishable from their eigenvalues A. Adding the "smoothing" term is expected 
not to alter the spatial localization of the filter. In practice, parameter a is 
selected to ensure 'computational uniqueness' of the smallest eigenvalue of 
D. Solutions obtained by the numerical implementation of the minimization 
of Cq are displayed in Figure U\ with various values for the smoothing pa- 
rameter a. Dashed lines correspond to the vector returned numerically as the 
"best" eigenvector of D (associated to the greatest eigenvalue), and the best 
eigenvector of D with parameter a chosen deliberately too small to ensure 
computationally uniqueness. Oscillating functions are indeed obtained. As a 
grows, the criterion selects non oscillating windows, two of which are shown 
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Fig. 7. Effect of the smoothing on the spectral and spatial shapes of PSWFs. 

by the plain lines. The loss measured by the increase of Cn is displayed in 
the legend of the lower panel and appears extremely small. In our example, 
the energy outside Q for the needlet built from C takes the value 2.78.10"^^, 
whereas its minimal possible value is 1.78.10"^^. 

A generalization of the Slepian concentration problem can be to consider other 
measures of concentration, such as L^, p = I, oo instead of L^. The criterion 
defined in Eq. (fT9l) becomes then 

*<^)--^ (21, 

where ||/||^ = /s |/(OI"de ifpG [l,oo) and ||/||- = esssup for a spher- 

ical function /. Unlike Slepian criterion Cq = Cq , these alternate criteria do 
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not lead to simple eigenvalue problems. They could be numerically optimized 
but this is beyond the scope of this paper. However we compare in Section [3?T] 
this criterion to the original one Cq. 



2. 3 Statistical criterion for optimal analysis with missing data 



Instead of focusing on the "geometrical" shape of the needlet, one may also 
optimize directly some alternate criterion of practical interest. 

In this section, we consider the following framework: given an underlying ran- 
dom field X on § to be analysed, a window function on § multiplying the 
field (for example a mask putting the field to zero in some regions) and a 
region D C § of interest in which the analysis is to be done, the aim is to 
get, in D, needlet coefficients of WX as close as possible to the coefficients 
computed from the uncorrupted field X. 

We shall assume statistical properties on the fields X, D and look for opti- 
mality of the filters on average. 

Assumption 2.1 (1) X is a real-valued Gaussian zero mean isotropic square 

integrable random field on S, with power spectrum (Ci). 
(2) W and D are deterministic elements o/H. 

Implicitly, X is a measurable mapping from some (A:',X, P) into (H, H), H 
being the Borel cx-filed of H. Let E denote the expectation operator under P. 
Recall that under Assumption 12.11 the covariance function on the field X is 
well defined and is given by 



It follows that E[X(0^] = (47r)^i EteN(2^ + ^)Ci- Moreover, the multipole 
moments (a^m) of X are complex Gaussian random variables. They are cen- 
tered, independent up to the relation aim = o.}-m satisfy Eda^oP) = 
iE(|a,^n = C,, m 7^ 0. 

Note that W and D can be indicator functions (binary masks) or any smooth 
functions on the sphere. 

A first attempt in this direction is the derivation of an unbiased estimate of 
the spectrum from the multipole moments and the empirical power spectrum 
of the weighted sky XW defined by 
2i+T respectively. It is well-known ( 



Peebles. 


1973; 


Hivon et al.. 


- 

2002 
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see also the compact proof in Appendix [A]) that 



2/" 4- 1 

E{Ce) = Mu'Ce with Mei' = E , (22) 

where the coefficients are defined by (IB.Sp . Note that the coupUng matrix 
M. depends on W only through its 'power spectrum' Cf^ . If M. is invertible, 
then {J^~^{Cii)) provides an unbiased estimate of (C^). 

Let now derive a criterion to design a window function h which minimises the 
effect of missing data in a needlet analysis procedure. We focus on a single band 
smoothed field {i.e. we fix one scale j) and the dependence on j is implicit 
in the notations. For a collection of couple of indices, say (^j, ^4)4=1 ... /, we 
use I](£i,mi),-i as a shorthand notation for the summation on ^i G N, mj G 

{ ^i) ■ ■ ■ ; ^j}; i 1, . . . , /. 

Given an analysis spectral window b = ih^^.^,--- ,be^^^) and its associated 
smoothing operator $ = Z]vin<^<^max ^^n^, the smoothed masked field is 

$xiy(o = E ^ J xiOwiOLei^ ■ ew. 

Write E[$X(0^] = {^tt)-^ Ee (^Ibj with a| = {2i + Let e denote the 

normalized diff^erence field 



Ei/2[$X(0 
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E ^eb',) E;,,™) bfa,raY,miO (23) 
where we have defined W = 1 — W, aim = {XW, V^m)- 

Suppose that {{{^k}k&Ki {^k)}keK) provides an exact Gauss quadrature for- 
mula at a degree 2£niax- Define Pk and the needlet coefficients of X and 
XW, respectively and define 

Pk — P'k 



Those random variables are normalized errors on the needlet coefficients in- 
duced by the application of the weight function W. If both X and XW are in 
H^maxi we easily check that ]E(/5|) = \/Xk{4:ny^ Y.ebjcrj and 

VA; eK, ek = e{^k) ■ 



The dispersion of either the continuous field e(^) or the finite set {ek}keK 
is taken as a measure of quality for an analysis $. This dispersion is not 
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measured on the whole sphere, since the difference e must be important in the 
regions where W is far from 1. In order to select the regions where e is to be 
minimized we introduce a function D = J2 demXim which provides a positive 
weight function in H. In the simplest case D can be Id for a region D of 
interest. More generally, D can be designed to give more or less importance to 
various regions of S according, for instance, to the need for reliability in the 
needlet coefficients. 



The coefficients or their continuous version e are used in two ways. The 
first one introduces a "tolerance" threshold a and counts the number of coef- 
ficients which are on average below this threshold. This m easure of the effi- 
ciency of a filter in t he presence of a mask is presented in iBaldi et al.l l2006l : 
Pietrobon et al.ll2006l but its optimization was not considered. The second one 
considers the integrated square error of e, weighted by the function D. It leads 
to a quadratic quadratic which is readily optimized. 



The first criterion, writes, for a binary function D, 

Efc:D(a.)=lP(|efc|<«) 



^{k : D{^k) = 1} 



(24) 



that is, the mean fraction of needlet coefficients corrupted by less than a 
normalized error a > 0. For an arbitrary function D, a possible generalization 
of ([24]) is 

^^^^^ ^ j:keKDi^k)n\ek\<a) 

In Subsection 13.21 we compare different windows using this criterion and a 
real mask. 

Alternately, consider now the mean integrated square error (MISE) 

i?(b)=E/D(0||e(Oirde (25) 

and define the optimal shape for the window b as 

b* = arg min R{h). (26) 

||b||=l 

Straightforward algebra leads to a close form expression of -R(b) depending 
on b, on the weight functions W and D, and on the power spectrum {C£)£(z^. 
Let wimydim denote the multipole coefficients of the weight functions W,D, 
respectively and 



m m' m" 



(see ( 1B.2[ ) for an expression as a function of the Wigner-3j coefficients). 
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Proposition 2.2 Under Assumption \2.1\ 



R{h) 



h'Qh 
h'ah 



where cr = diag((cr|)) and Q is the matrix with entries 

^ — ' ^ — '(ti,rni) ^ 
m,m' 

^ ,m,,)..2,3,4 "^^^"^2 "^^^ams de^m. 



il i2 i 






* 




mi m2 m 




mi rris m' 




m 7714 fn' 



■ (27) 



// both W and D are axisymmetric the ten-tuple summations above reduce to 
a five-tuple one 



Qw — X! X! Cij^Wi^QWi^odi^fi 

m £1/2/3/4 













£1 £2 £ 




£1 £3 £' 




£ £4 £' 


771 771 




m m 




m m 



Wm 



In the next section we shall give some illustrative examples of optimal spectral 
windows h* in the particular axisymmetric case. 

Remark 2.3 As in the Slepian's problem, the design of an optimal filter re- 
duces to an eigenvalue problem. In particular, if > for any £ G L, write 
h\ = a (hi. Let be an eigenvector associated with the lowest eigenvalue of 
Q\ q\^, = {(Ti(Te)-^Qu'- Then b* := crbt7||crbt*|| is a solution of 

Remark 2.4 For those sums to be tractable, one has to assume that D, W , 
Ce have finite support in the frequency domain, i.e. that the windows D and 
W are smooth (or apodized) and = for large enough £. 

Remark 2.5 The matrix Q being a second-order moment for the random field 
X, it can also be approximated by a moment estimator using Monte-Carlo 
experiments. This remark is of important practical interest as we are mostly 
concerned with non zonal masks. 
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3 Examples, numerical results 



3. 1 Comparison of filters for various criteria 



In Section [21 we considered several criteria measuring the localization prop- 
erties of filters, and derived explicit or computational optimization for some 
of them. In Table [U we compare the scores reached by the filters displayed 
in Figure [H The columns indexed by list the values Cq^^ (ip) defined in 
Eq. (fT9l) . More generally, the columns indexed by 17-9 correspond to the val- 
ues Cf^gi'ip) defined in Eq. ( |2T1) . A column lists the values of 1 — E{a) defined 
in Eq. (i24ll . applied with the mask KpO of Fig. [8] and a tolerance parameter 
a = 10% (see next subsection for more details). A last column gives, by way of 
illustration only, the value of the "uncertainty product" A^{ip) x Al(^), where 



— II f Ei>(E)^d^P 



Narcowich and WardI ( 1996h proved that A^{ip) x Al(^) > 1. 



The PSWFs perform the best not only for the criterion which they optimize, 
but also in most cases for the criteria where the norm is replaced by ones, 
p = 1 and p = oo, with the same opening angles Oq. Although the KpO mask 
has many small cut areas all over the sphere, most of the 11 filters presented 
here allow to retain more than 60% of the outside-mask coefficients Pk if a 
10% error due to the presence of the mask is accepted. The performance w.r.t. 
this criterion goes up to 85% for the PSWF optimally concentrated in a cap 
of 1°. However, the choice of arbitrary value of a has a major impact on the 
ranking of the filters. This point is investigated in the next subsection. 



3.2 Robustness of needlets coefficients 



In this Subsection, we illustrate the performances of various window functions 
using the criterion (l24l) . We have run = 30 Monte-Carlo experiments to 
estimate the numerator of Ey^{a). The random fields X are drawn using the 
(C/l-s pectrum of the best-fit ting model for the CMB estimated by the WMAP 



team ( IHinshaw et al.l . l2006l ). The mask W was chosen as KpO, displayed in 
Figure [HI which masks the galactic plane and many point sources. The band 
is L = [256,1024]. 

Figure [9] compares the increasing functions -E'b(') corresponding to various 
filters b. There is no "uniformly best" {i.e. highest in the figure) needlet: some 
allow to retain more coefficients when the constraint imposed on the error is 
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l.Oe-01 


1.5e-02 


5.0e-02 


1.9e-02 


5.1e-03 


6.4e-05 


2.6e-01 


2.7 


Spline, order 7 


4.0e-02 


1.3e-02 


2.0e-03 


4.8e-08 


5.0e-01 


2.9e-01 


1.3e-01 


1.7e-03 


6.0e-02 


2.7e-02 


7.1e-03 


1.2e-05 


3.3e-01 


3.1 


Spline, order 15 


6.1e-02 


2.5e-02 


4.9e-03 


4.0e-07 


5.9e-01 


4.0e-01 


2.2e-01 


2.3e-03 


6.9e-02 


3.3e-02 


9.8e-03 


7.0e-05 


4.1e-01 


3.7 


Spline, order 21 


7.2e-02 


3.1e-02 


7.1e-03 


7.7e-06 


6.2e-01 


4.5e-01 


2.7e-01 


l.Oe-02 


7.3e-02 


3.7e-02 


l.le-02 


2.7e-04 


4.6e-01 


4.1 


Prolate, cap 0.5° 


1.2e-02 


6.0e-03 


3.4e-03 


9.5e-04 


8.5e-01 


8.2e-01 


8.0e-01 


7.2e-01 


5.1e-02 


l.Oe-02 


5.6e-03 


l.Oe-03 


6.5e-01 


9.8 


Prolate, cap 1° 


6.7e-02 


4.3e-05 


5.8e-06 


1.7e-06 


3.8e-01 


l.,3e-01 


1.2e-01 


l.le-01 


l.le-01 


2.0e-03 


2.0e-04 


5.0e-05 


1.5e-01 


3.1 


Prolate, cap 1.5° 


1.2e-01 


1.5e-03 


3.4e-07 


1.2e-08 


4.3e-01 


5.3e-02 


l.Oe-02 


8.8e-03 


1.3e-01 


1.7e-02 


1.4e-04 


4.5e-06 


1.7e-01 


3.6 


Prolate, cap 5° 


l.le-01 


6.7e-03 


6.5e-04 


5.7e-14 


5.0e-01 


1.8e-01 


6.8e-02 


l.le-06 


1.2e-01 


2.2e-02 


5.9e-03 


2.6e-08 


2.4e-01 


3.6 


Exponential 


1.8e-02 


3.2e-03 


l.Oe-03 


l.Oe-05 


4.4e-01 


2.6e-01 


1.9e-01 


4.8e-02 


4.4e-02 


1.4e-02 


5.7e-03 


1.9e-04 


2.7e-01 


2.7 


B-Spline 


l.le-02 


1.3e-03 


3.9e-04 


1.3e-05 


4.8e-01 


3.3e-01 


2.7e-01 


1.5e-01 


3.1e-02 


6.8e-03 


2.5e-03 


1.5e-04 


2.1e-01 


1.2 


Mexican hat 


6.4e-01 


l.le-02 


8.5e-07 


7.,3e-12 


7.9e-01 


8.8e-02 


7.1e-04 


1.6e-04 


4.7e-01 


7.9e-02 


4.5e-04 


1.4e-07 


4.9e-01 


3.0 



Table 1. Comparison of the eleven filters of Fig. HJ the nine first of which are band-limited in L=[256,1024]. 



Fig. 8. KpO mask. 



loose enough, but their efficiency decreases faster as a goes to zero. Inspect 
e.g. the PSWF family. 




Fig. 9. Proportion -E'b(a) of coefficients uncontaminated at tolerance level a. 



3.3 Some MISE-optimal filters for axisymmetric weight functions 



We present here the results of the optimization fl26l) in the case of axisymmetric 
weight functions W. For simplicity, the reconstruction weight function D is 
taken equal to W. We stick to the CMB spectrum of previous subsection. 

Figure [TOl displays some of the masks W used in the experiments. The apodiza- 
tion in simply a cosine- arch junction between and 1, on a 2-degrees angular 
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(b) 




(c) 



(d) 



Fig. 10. Four different apodized masks. The degree of apodization, measured as the 
width of the cosine-arch 0-1 junction is, 2 degrees. 



range. This means that the data is available on the dark regions, and that its 
L = [imin, ^max]-band-limited part has to be recovered in this area too. 



On Figure \TT\ we have plotted the optimal filter in the i?(b)-sense for the 
masks of Figure fTOl together with different PSWFs. The criterion captures the 
symmetry of the mask (a) (the shape of the matrix Q is a "checkerboard"), 
and the optimal filter is thus zero on all even (here) or all odd multipoles. The 
associated axisymmetric needlet ifj is symmetric w.r.t. the equatorial plane, 
and thus is well concentrated around both the North and the South poles. 
Such solutions are very sensitive to the modifications of the masks. 



We conducted a small Monte-Carlo study to confirm the benefit of our ap- 
proach. We have compared our best filters b* to PSWFs with different open- 
ing. On Figure [121 we show the box-plots of the distribution of the statistic 
i?(h) for all those filters. Stars are plotted at the position of the estimated 
value of Ei?(b) and the horizontal line is this value for b*. The right vertical 
scale is for the relative error (in percent) with respect to Ei?(b*). 



Fig. 12(a) illustrates the strong benefit of a filter that captures the geometry of 
the mask. The relative improvement with respect to the best PSWF is of order 
20%. It should be noted however that the shape of this optimal filter (described 
above) may lead to a misleading space-frequency picture. In some other cases, 
as shown in Figure 12(b), the relative improvement from the best PSWF to 
the best filter at all is very slight (a few percents). Here, the most favorable 
feature of our approach is that there is no tuning parameters (opening of the 
PSWF for instance, or the order of the splines window functions if they are 
taken as alternatives) to be found before the analysis. 
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Multipole 



Fig. 11. Shape of optimal window functions (plain lines) and PSWF (coloured and 
dashed lines) with various openings. 
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Fig. 12. Mean square error in analysis. Stars are potted at the estimated values for 
R(h). For Figure 12(a) , the mask is Figure [TOlfa) and L 
the mask is Figure [TO] and L = [20, 30] 



[5,15]. For Figure 12(b) 



25 



4 Conclusion 



A flexible way of analysing a fleld on the sphere in a space-f requency manner 



has be en presented. It is based on the needlet construction of lNarcowich et al. 



( 120061 ). The proposed analysis functions form a frame in the space a square- 
integrable functions on the sphere. Decompositions are essentially operating 
in the Spherical Harmonics domain, leading to fast implementations. Various 
criteria are used to design good spectral windows. This optimization can lead 
to decisive improvement in high precision applications such as modern cosmol- 
ogy (CMB spectral estimation, component separation, etc.), where localized 
analysis is crucial. 
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A Proofs 



Proof of Propositions 11.41 and 11.51 Propositions 11.41 is a particular case 
of Proposition O Indeed ©-jHl) imply ^ with Ci = C2 = 1- Together 
with ffTOl) we get ^l^^ = (3^'' and ■i/'^"'^ = ip\^\ Prove now Proposition 11.51 
Firstly, using successively ^ and the quadrature assumption (remind that 
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for any £,f < d, (n,X)(n,X) G ®ft,m, 



E m' = j: E 



= E E A? 

= E E ^^^^ E A?n,x(a)n,,x(a) 

= E E / n,x(e)n,,x(Ode 

= E E ^^^^s^i' / in.x(0|Me 
= EE#)'lin.xf . 

Using m and = E^lin.Xlp, we get Ci||Xf < < C2\\Xf. 

Prove now that {ipj^^) is the dual frame of (ipk^). Write 

e'=oe=o ^ 

^=0 



Then, for any j eN, ke K^^^ 



j',k' e=o e'=o 



(A?)^/^ E E E fe/^ E ^''L^iiP ■ ^')L.{^' 
(Ai^'V/^ E E ^^"^4^'^ E L Le{^' ■ ■ Odr 



j'ej £=0 



£'=0 



£=0 j'ej 
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The asserti o ns (fT3l) are a consequence of the dual frame property (see e.g. 
Daubechiesl . Il992l ). 



Proof of Proposition 11.61 From Definition 11.31 of ipk and Eq ([T]) 



□ 



Proof of Eq. (El) 

(2£ + l)E(Q) = 



£ 

E // E{x{ox{e)}Y,uoyim{e)w{ow{ad^de 

ff {E Ce^L.i^ ■ ■ OWiOWiOd^de 

II L,,{i-OW{0W{Od^di' 



£'gn o<i"<e+e' 

E a«'£"(2r + l)C,'^ 

£'GN 0<l"<l+l' 



Proof of Proposition[2l2] As XW = E(fi,mi) a^imiy^imi ^(£2,^2) W£2m2^^2m2, 

a^m := {XW,Yim) = E^^ 1 «^imiW^2m2 

V-f-z — 1,2 

Together with E[a£ma|/^/] = Ci5ui5mm' it yields 



mi m2 m 



¥\ai>mO*ern'] " E(£,,mi)i=i 2 3 *^^l^^2m2 ^^3^3 



£1 4 ^ 




^1 4 e 


mi m2 m 




mi m3 m' 



(A.l) 



Combining ( !23l) and (lA.ll ) we get 



^ E E ^K"i«Iw]E(£,^,)^^4m4 / ^^4m4(0^^m(0^/m'(^^ 

£GN £/'GN mm' ^ ' 



ij^^e^l) ^ E hpheQw- 

teN £,£'gn 
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If W is axisymmetric, 



^2 A 



where we used the fact that 



nil m 



iTLi m' 









£l i2 i 






m m 




m m 



Ml'ra (A. 2) 



if m 7^ m" and that tfiyo and 



m m 



m m" 

are real. If D is axisymmetric and with -D«'m '■= ^^£4 ^^^4,0 



m m 



(4,mi)«=i,2,3 
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□ 



B Legendre polynomials, spherical harmonics and related useful 
formulae 



Usually, Pi{z) denotes the Legendre polynomial of order £, normalized by 
Pi{l) = 1. For our purposes, it is more convenient to use a different normal- 
ization 



LAz) 



Air 



because we get coefficient-free properties like 



and 



m=-i 



Js 



(B.l) 



In other words, Li is the polynomial kernel of the harmonic projection on H^. 
We have f^^ Pi{zfdz = and J^^ Lt{zfdz - ^ 



Stt^ • 



The spherical harmonics are explicitly given in a factorized form in terms of 
the associated Legendre polynomials and the complex exponentials as 



(2f + 1) {e-m)\ 
^ Air {i + m)\ 



Pimicos6)e 



irrup 
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where Pemix 



The following equations relate the integral of the product of three complex 
spherical harmonics over the total solid angle or three Legendre polynomials 
with the Wign er-3? coefficients (for a d efinition in terms of Clebsh-Gordan 
coefficients, see IVarshalovich et alJ (If 9881 ). pp235-). 



12 t3 

mi m2 rris 



-1 



/(2£i + l)(2f2 + l)(2£3 + l) i h ^2 h 



I 
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-1 ^2 



ymi 7712 —m^ 
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Le{z)U>{z')Le'iz")dzdz'dz' 



,^ „ (2£ + l)(2f + l)(2r + l) 



(4vr) 



(47r)= 



t2 t3 





From (IB.3I) and / L^L^/ = 5^/'^^ we get: 



0<i"<£+t' 



with 



(2£+l)(2f + 1) 
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(B.4) 



(B.5) 
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